package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.*;

import java.util.Arrays;
import java.util.HashSet;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.*;

public class FES2Q41 extends FES {

    public FES2Q41(Mesh mesh2D) {
        super(mesh2D, 1);
        eleStiff = FES2Q41::elementStiff;
        eleMass = FES2Q41::elementMass;
        eleSource = FES2Q41::elementSource;
        eleFlux = FES2T31::elementFlux;
    }

    public double getValue(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
        return getValue0(u, xyz, eleIdx, eleVal, CG::pointInQ4);
    }

    ;

    public double getValue(double[] u, double[] xyz) {
        return getValue0(u, xyz, CG::pointInQ4);
    }

//	public Pair<Double, Integer> getValue(double[] u, double[] xy, int previousEleIdx, double[][] coef, int[] eleIdx, double[] eleVal) {
//		ParamCheck.checkEqual(u.length, this.ndof);
//		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
//		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
//		double value = 0;
//		double x = xy[0], y = xy[1];
//		double[] v = new double[] {1.0, x, y, x * y};
//		double[][] coord = new double[4][2];
//		int[] ele = mesh.getElement(previousEleIdx);
//		mesh.getCoordPerElement(ele, coord);
//		if (CG.pointInQ4(coord, xy)) {
//			
//			MVO.matmul(coef, v, eleVal);
//			//CG.Q4ShapeFunction(coord, xy, eleVal);
//			for (int j = 0; j < 4; j++) {
//				//eleIdx[j] = ele.idx[j];
//				value += u[ele[j]] * eleVal[j];
//			}
//			return new Pair<Double, Integer>(value, previousEleIdx);
//		} else {
//			for (int i = 0, n = mesh.getNt(); i < n; i++) {
//				ele = mesh.getElement(i);
//				mesh.getCoordPerElement(ele, coord);
//				if (CG.pointInQ4(coord, xy)) {
//					//CG.Q4ShapeFunction(coord, xy, eleVal);
//					CG.Q4ShapeFuncionCoef(coord, coef);
//					MVO.matmul(coef, v, eleVal);
//					for (int j = 0; j < 4; j++) {
//						eleIdx[j] = ele[j];
//						value += u[ele[j]] * eleVal[j];
//					}
//					return new Pair<Double, Integer>(value, i);
//				}
//			}
//			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
//		}
//
//	}

    public void getValue(double[] u, double[] xy, DoubleIntegerPair data, double[][] coef, int[] eleIdx,
                         double[] eleVal) {
        ParamCheck.checkEqual(u.length, this.ndof);
        ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
        ParamCheck.checkEqual(eleVal.length, this.nPerEle);
        double value = 0;
        double x = xy[0], y = xy[1];
        double[] v = new double[]{1.0, x, y, x * y};
        double[][] coord = new double[4][2];
        int[] ele = mesh.getElement(data.second);
        mesh.getCoordPerElement(ele, coord);
        if (CG.pointInQ4(coord, xy)) {

            MVO.matmul(coef, v, eleVal);
            // CG.Q4ShapeFunction(coord, xy, eleVal);
            for (int j = 0; j < 4; j++) {
                // eleIdx[j] = ele.idx[j];
                value += u[ele[j]] * eleVal[j];
            }
            data.first = value;
            return;
        } else {
            for (int i = 0, n = mesh.getNt(); i < n; i++) {
                ele = mesh.getElement(i);
                mesh.getCoordPerElement(ele, coord);
                if (CG.pointInQ4(coord, xy)) {
                    // CG.Q4ShapeFunction(coord, xy, eleVal);
                    CG.Q4ShapeFuncionCoef(coord, coef);
                    MVO.matmul(coef, v, eleVal);
                    for (int j = 0; j < 4; j++) {
                        eleIdx[j] = ele[j];
                        value += u[ele[j]] * eleVal[j];
                    }
                    // data.setFirst(value);
                    // data.setSecond(i);
                    data.first = value;
                    data.second = i;
                    return;
                }
            }
            throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
        }
    }

    public void getValue_Square(double[] u, double[] xy, DoubleIntegerPair data, double[][] coef, int[] eleIdx,
                                double[] eleVal) {
        ParamCheck.checkEqual(u.length, this.ndof);
        ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
        ParamCheck.checkEqual(eleVal.length, this.nPerEle);
        double value = 0;
        double x = xy[0], y = xy[1];
        double[] v = new double[]{1.0, x, y, x * y};
        double[][] coord = new double[4][2];
        int[] ele = mesh.getElement(data.second);
        mesh.getCoordPerElement(ele, coord);
        if (CG.pointInSquare(coord, xy)) {

            MVO.matmul(coef, v, eleVal);
            // CG.Q4ShapeFunction(coord, xy, eleVal);
            for (int j = 0; j < 4; j++) {
                // eleIdx[j] = ele.idx[j];
                value += u[ele[j]] * eleVal[j];
            }
            data.first = value;
            return;
        } else {
            for (int i = 0, n = mesh.getNt(); i < n; i++) {
                ele = mesh.getElement(i);
                mesh.getCoordPerElement(ele, coord);
                if (CG.pointInSquare(coord, xy)) {
                    // CG.Q4ShapeFunction(coord, xy, eleVal);
                    CG.Q4ShapeFuncionCoef(coord, coef);
                    MVO.matmul(coef, v, eleVal);
                    for (int j = 0; j < 4; j++) {
                        eleIdx[j] = ele[j];
                        value += u[ele[j]] * eleVal[j];
                    }
                    // data.setFirst(value);
                    // data.setSecond(i);
                    data.first = value;
                    data.second = i;
                    return;
                }
            }
            throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
        }
    }

    /**
     * 计算单元热刚度矩阵
     *
     * @param xy        单元结点坐标数组， 行数为nPerElement, 列数为2
     * @param kxx       单元热传导系数, 行数为nPerElement, 列数为2
     * @param tp        热传导问题，可取为COMMON, 与AXISYMMETRIC
     * @param eleMatrix 存储计算得到的单元热刚度矩阵
     */
    public static void elementStiff(int[] ele, double[][] xy, double[][] kxx, BVPType tp, double[][] eleMatrix) {
        double[] x = MVO.col(xy, 0);
        double[] y = MVO.col(xy, 1);
        double[] dx = new double[4];
        double[] dy = new double[4];
        double[] v = new double[4];
        double dxU, dyU;
        double c11, c22;
        double[][] temp;
        double[] N = new double[3];
        double[] coef0_11, coefX_11, coefX2_11, coefX3_11, coefU_11, coefU2_11, coefU3_11;
        double[] coef0_22, coefX_22, coefX2_22, coefX3_22, coefU_22, coefU2_22, coefU3_22;
        double[] coefdX_11, coefdX2_11, coefdX3_11, coefdY_11, coefdY2_11, coefdY3_11;
        double[] coefdX_22, coefdX2_22, coefdX3_22, coefdY_22, coefdY2_22, coefdY3_22;
        double[] u;
        double weight;
        double xi, et, detJ = 0.0, wxarea = 0.0, c;
        MVO.fill(eleMatrix, 0.0);
        switch (tp) {
            case COMMON:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(kxx[0], v);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += (dx[si] * dx[sj] + dy[si] * dy[sj]) * wxarea;
                            }
                        }
                    }
                }
                // System.out.println(wxarea + "\n" + MVO.toString(eleMatrix));
                break;
            case AXISYMMETRIC:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(kxx[0], v) * MVO.dot_product(x, v);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += (dx[si] * dx[sj] + dy[si] * dy[sj]) * wxarea;
                            }
                        }
                    }
                }
                break;
            case NONLINEAR_POLYNORMIAL_2:
                //a11 + f11 u + g11 u^2 + h11 u^3
                //a22 + f22 u + g22 u^2 + h22 u^3
                coef0_11 = kxx[0];
                coefU_11 = kxx[1];
                coefU2_11 = kxx[2];
                coefU3_11 = kxx[3];
                coef0_22 = kxx[4];
                coefU_22 = kxx[5];
                coefU2_22 = kxx[6];
                coefU3_22 = kxx[7];
                u = kxx[8];
                dxU = MVO.dot_product(u, dx);
                dyU = MVO.dot_product(u, dy);
                for (int i = 0; i < 3; i++) {
                    xi = GP3[i];
                    for (int j = 0; j < 3; j++) {
                        et = GP3[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        double meanU = MVO.dot_product(u, v);
                        c11 = MVO.dot_product(coefU_11, v);
                        c11 += 2 * meanU * MVO.dot_product(coefU2_11, v);
                        c11 += 3 * meanU * meanU * MVO.dot_product(coefU3_11, v);
                        c22 = MVO.dot_product(coefU_22, v);
                        c22 += 2 * meanU * MVO.dot_product(coefU2_22, v);
                        c22 += 3 * meanU * meanU * MVO.dot_product(coefU3_22, v);
                        wxarea = GW3[i] * GW3[j] * detJ;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += (c11 * dxU * dx[si] * v[sj] + c22 * dyU * dy[si] * v[sj]) * wxarea;
                            }
                        }
                    }
                }
            case NONLINEAR_POLYNORMIAL_1:
                //a11 + f11 u + g11 u^2 + h11 u^3
                //a22 + f22 u + g22 u^2 + h22 u^3
                coef0_11 = kxx[0];
                coefU_11 = kxx[1];
                coefU2_11 = kxx[2];
                coefU3_11 = kxx[3];
                coef0_22 = kxx[4];
                coefU_22 = kxx[5];
                coefU2_22 = kxx[6];
                coefU3_22 = kxx[7];
                u = kxx[8];
                for (int i = 0; i < 3; i++) {
                    xi = GP3[i];
                    for (int j = 0; j < 3; j++) {
                        et = GP3[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c11 = MVO.dot_product(coef0_11, v);
                        c22 = MVO.dot_product(coef0_22, v);
                        double meanU = MVO.dot_product(u, v);
                        c11 += MVO.dot_product(coefU_11, v) * meanU;
                        c22 += MVO.dot_product(coefU_22, v) * meanU;
                        c11 += MVO.dot_product(coefU2_11, v) * meanU * meanU;
                        c22 += MVO.dot_product(coefU2_22, v) * meanU * meanU;
                        c11 += MVO.dot_product(coefU3_11, v) * meanU * meanU * meanU;
                        c22 += MVO.dot_product(coefU3_22, v) * meanU * meanU * meanU;
                        wxarea = GW3[i] * GW3[j] * detJ;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += (c11 * dx[si] * dx[sj] + c22 * dy[si] * dy[sj]) * wxarea;
                            }
                        }
                    }
                }
                break;
            case NONLINEAR_DERIVATIVE_2:
                //a11 + b11 du/dx + c11 (du/dx)^2 + d11 (du/dx)^3 + f11 du/dy + g11 (du/dy)^2 + h11 (du/dy)^3
                //a22 + b22 du/dx + d22 (du/dx)^2 + d22 (du/dx)^3 + f22 du/dy + g22 (du/dy)^2 + h22 (du/dy)^3
                double c12, c21;
                coef0_11 = kxx[0];
                coefdX_11 = kxx[1];
                coefdX2_11 = kxx[2];
                coefdX3_11 = kxx[3];
                coefdY_11 = kxx[4];
                coefdY2_11 = kxx[5];
                coefdY3_11 = kxx[6];
                coef0_22 = kxx[7];
                coefdX_22 = kxx[8];
                coefdX2_22 = kxx[9];
                coefdX3_22 = kxx[10];
                coefdY_22 = kxx[11];
                coefdY2_22 = kxx[12];
                coefdY3_22 = kxx[13];
                u = kxx[14];
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        dxU = MVO.dot_product(dx, u);
                        dyU = MVO.dot_product(dy, u);
                        c11 = (MVO.dot_product(coefdX_11, v) + (2 * MVO.dot_product(coefdX2_11, v) + 3 * MVO.dot_product(coefdX3_11, v) * dxU) * dxU) * dxU;
                        c12 = (MVO.dot_product(coefdY_11, v) + (2 * MVO.dot_product(coefdY2_11, v) + 3 * MVO.dot_product(coefdY3_11, v) * dyU) * dyU) * dxU;
                        c21 = (MVO.dot_product(coefdX_22, v) + (2 * MVO.dot_product(coefdX2_22, v) + 3 * MVO.dot_product(coefdX3_22, v) * dxU) * dxU) * dyU;
                        c22 = (MVO.dot_product(coefdY_22, v) + (2 * MVO.dot_product(coefdY2_22, v) + 3 * MVO.dot_product(coefdY3_22, v) * dyU) * dyU) * dyU;
                        wxarea = GW2[i] * GW2[j] * detJ;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += (c11 * dx[si] * dx[sj] + c12 * dx[si] * dy[sj]
                                        +c21 * dy[si] * dx[sj] + c22 * dy[si] * dy[sj]) * wxarea;
                            }
                        }
                    }
                }
            case NONLINEAR_DERIVATIVE_1:
                coef0_11   = kxx[0];
                coefdX_11  = kxx[1];
                coefdX2_11 = kxx[2];
                coefdX3_11 = kxx[3];
                coefdY_11  = kxx[4];
                coefdY2_11 = kxx[5];
                coefdY3_11 = kxx[6];
                coef0_22   = kxx[7];
                coefdX_22  = kxx[8];
                coefdX2_22 = kxx[9];
                coefdX3_22 = kxx[10];
                coefdY_22  = kxx[11];
                coefdY2_22 = kxx[12];
                coefdY3_22 = kxx[13];
                u = kxx[14];
                for (int i = 0; i < 3; i++) {
                    xi = GP3[i];
                    for (int j = 0; j < 3; j++) {
                        et = GP3[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        dxU = MVO.dot_product(u, dx);
                        dyU = MVO.dot_product(u, dy);
                        c11 = MVO.dot_product(coef0_11, v);
                        c11 += MVO.dot_product(coefdX_11, v)  * dxU;
                        c11 += MVO.dot_product(coefdX2_11, v) * dxU * dxU;
                        c11 += MVO.dot_product(coefdX3_11, v) * dxU * dxU * dxU;
                        c11 += MVO.dot_product(coefdY_11, v)  * dyU;
                        c11 += MVO.dot_product(coefdY2_11, v) * dyU * dyU;
                        c11 += MVO.dot_product(coefdY3_11, v) * dyU * dyU * dyU;
                        c22 = MVO.dot_product(coef0_22, v);
                        c22 += MVO.dot_product(coefdX_22, v)  * dxU;
                        c22 += MVO.dot_product(coefdX2_22, v) * dxU * dxU;
                        c22 += MVO.dot_product(coefdX3_22, v) * dxU * dxU * dxU;
                        c22 += MVO.dot_product(coefdY_22, v)  * dyU;
                        c22 += MVO.dot_product(coefdY2_22, v) * dyU * dyU;
                        c22 += MVO.dot_product(coefdY3_22, v) * dyU * dyU * dyU;
                        wxarea = GW3[i] * GW3[j] * detJ;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += (c11 * dx[si] * dx[sj] + c22 * dy[si] * dy[sj]) * wxarea;
                            }
                        }
                    }
                }

                break;

            default:
                throw new IllegalArgumentException("Illegal BVPType!");
        }

        // System.out.println(MVO.toString(eleMatrix));
        // System.exit(0);
    }

    public static void elementMass(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleMatrix) {
        double[] x = MVO.col(xy, 0);
        double[] y = MVO.col(xy, 1);
        double[] dx = new double[4];
        double[] dy = new double[4];
        double[] v = new double[4];
        double xi, et, detJ, wxarea, c;
        MVO.fill(eleMatrix, 0.0);
        switch (tp) {
            case CONSIST_MASS:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(coef[0], v);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += v[si] * v[sj] * wxarea;
                            }
                        }
                    }
                }
                break;
            case LUMP_MASS:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(coef[0], v);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            eleMatrix[si][si] += v[si] * wxarea;
                        }
                    }
                }
                break;
            case AXIS_CONSIST_MASS:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(coef[0], v) * MVO.dot_product(x, v);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            for (int sj = 0; sj < 4; sj++) {
                                eleMatrix[si][sj] += v[si] * v[sj] * wxarea;
                            }
                        }
                    }
                }
                break;
            case AXIS_LUMP_MASS:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(coef[0], v) * MVO.dot_product(x, v);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            eleMatrix[si][si] += v[si] * wxarea;
                        }
                    }
                }
                break;
            default:
                throw new IllegalArgumentException("Illegal BVPType!");
        }
    }

    public void assembleHeatThetaScheme(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix,
                                        SRMatrix RMatrix) {
        double[][] kxx = new double[1][nPerEle];
        double[][] rho = new double[1][nPerEle];
        double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
        double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
        double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
        double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
        double dt = param[0];
        double theta = param[1];
        for (int i = 0; i < nE; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            int label = mesh.getElementLabel(i);
            for (int j = 0, nj = nPerEle; j < nj; j++) {
                kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
                rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
            }
            elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
            elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
            for (int ii = 0; ii < dofPerElement; ii++) {
                for (int jj = 0; jj < dofPerElement; jj++) {
                    eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
                    eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
                }
            }
            LMatrix.assemble(ele, eleLMatrix);
            RMatrix.assemble(ele, eleRMatrix);
        }
    }

    public void assembleHeatThetaScheme(double[] coefFunc, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
        double[][] kxx = new double[1][nPerEle];
        double[][] rho = new double[1][nPerEle];
        double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
        double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
        double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
        double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
        Arrays.fill(kxx[0], coefFunc[0]);
        Arrays.fill(rho[0], coefFunc[1]);
        double dt = coefFunc[2];
        double theta = coefFunc[3];
        for (int i = 0; i < nE; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
            elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
//			System.out.println(MVO.toString(eleMatrixMass));
//			System.out.println(MVO.toString(eleMatrixStiff));
//			System.exit(0);
            for (int ii = 0; ii < dofPerElement; ii++) {
                for (int jj = 0; jj < dofPerElement; jj++) {
                    eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
                    eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
                }
            }
            LMatrix.assemble(ele, eleLMatrix);
            RMatrix.assemble(ele, eleRMatrix);
        }
    }

    public static void elementSource(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
        double[] x = new double[]{xy[0][0], xy[1][0], xy[2][0], xy[3][0]};
        double[] y = new double[]{xy[0][1], xy[1][1], xy[2][1], xy[3][1]};
        double[] dx = new double[4];
        double[] dy = new double[4];
        double[] v = new double[4];
        double xi, et, detJ, wxarea, c;
        MVO.setValue(eleF, 0.0);
        switch (tp) {
            case COMMON:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(f[0], v);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            eleF[si] += v[si] * wxarea;
                        }
                    }
                }
                break;
            case AXISYMMETRIC:
                for (int i = 0; i < 2; i++) {
                    xi = GP2[i];
                    for (int j = 0; j < 2; j++) {
                        et = GP2[j];
                        detJ = CG.ISOMap2D(CG::Q4ShapeFunction, x, y, xi, et, v, dx, dy);
                        c = MVO.dot_product(f[0], v) * MVO.dot_product(v, x);
                        wxarea = GW2[i] * GW2[j] * detJ * c;
                        for (int si = 0; si < 4; si++) {
                            eleF[si] += v[si] * wxarea;
                        }
                    }
                }
                break;
            default:
                throw new IllegalArgumentException("Illegal BVPType!");
        }
    }

    public void assembleFlux(ScalarFunc[] coeffFunc, double[] param, BVPType tp, double[] RHS, int... bdLabel) {
        assembleBoundaryVector(coeffFunc, param, FES2T31::elementFlux, tp, RHS, bdLabel);
    }

    public void assembleFlux(double[] constCoef, BVPType tp, double[] RHS, int... bdLabel) {
        assembleBoundaryVector(constCoef, FES2T31::elementFlux, tp, RHS, bdLabel);
    }

    public void assembleFlux(double[][] vecCoef, BVPType tp, double[] RHS, int... bdLabel) {
        assembleBoundaryVector(vecCoef, FES2T31::elementFlux, tp, RHS, bdLabel);
    }

    public void assembleBoundaryMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A,
                                     int... bdLabel) {
        assembleBoundaryMatrix(coefFunc, param, FES2T31::elementBoundaryMass, tp, A, bdLabel);
    }

    public void assembleBoundaryMass(double[] constConvect, BVPType tp, SRMatrix A, int... bdLabel) {
        assembleBoundaryMatrix(constConvect, FES2T31::elementBoundaryMass, tp, A, bdLabel);
    }

    public void assembleBoundaryMass(double[][] vectorConvect, BVPType tp, SRMatrix A, int... bdLabel) {
        assembleBoundaryMatrix(vectorConvect, FES2T31::elementBoundaryMass, tp, A, bdLabel);
    }

    public void computeGradient(double[] u, double[][] gradU) {
        double xi, et;
        double[] uu = new double[4];
        double[] v = new double[4];
        double[] dx = new double[4];
        double[] dy = new double[4];
        double[] count = new double[u.length];
        for (int k = 0, n = mesh.getNt(); k < n; k++) {
            int[] ele = mesh.getElement(k);
            mesh.getCoordPerElement(ele, coordRegion);
            for (int i = 0; i < 4; i++) {
                uu[i] = u[ele[i]];
            }
            double[] dxUU = new double[4];
            double[] dyUU = new double[4];
            double detJ = 0.0;
            for (int i = 0; i < 2; i++) {
                xi = GP2[i];
                for (int j = 0; j < 2; j++) {
                    et = GP2[j];
                    detJ = CG.ISOMap2D(CG::Q4ShapeFunction, MVO.col(coordRegion, 0), MVO.col(coordRegion, 1), xi, et, v,
                            dx, dy);
                    dxUU[j + i * 2] = MVO.dot_product(uu, dx) * detJ;
                    dyUU[j + i * 2] = MVO.dot_product(uu, dy) * detJ;
                }
            }
            double[] temp = MVO.matmul(Constants.GaussMatrix2D, dxUU);
            double[] temp2 = MVO.matmul(Constants.GaussMatrix2D, dyUU);
            for (int i = 0; i < 4; i++) {
                gradU[0][ele[i]] += temp[i];
                gradU[1][ele[i]] += temp2[i];
                count[ele[i]] += detJ;
            }
        }

        for (int i = 0; i < mesh.getNv(); i++) {
            gradU[0][i] /= count[i];
            gradU[1][i] /= count[i];
        }
    }

    public void computeGradient2(double[] u, double[][] gradU) {
        int len = u.length;
        ParamCheck.checkEqual(len, this.ndof);
        ParamCheck.checkEqual(len, gradU[0].length);
        ParamCheck.checkEqual(len, gradU[1].length);

        double[] angle = new double[4];
        double[][] coord = new double[4][2];
        double[][] grad = new double[4][2];
        double[] angleCount = new double[this.GetNdof()];
        double[] eleU = new double[4];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coord);
            eleU[0] = u[ele[0]];
            eleU[1] = u[ele[1]];
            eleU[2] = u[ele[2]];
            eleU[3] = u[ele[3]];
            CG.gradientValue2DQ4(coord, eleU, grad);
            CG.angle2DQ4(coord, angle);
            for (int j = 0; j < 4; j++) {
                gradU[0][ele[j]] += grad[j][0] * angle[j];
                gradU[1][ele[j]] += grad[j][1] * angle[j];
                angleCount[ele[j]] += angle[j];
            }
        }
        for (int i = 0, n = mesh.getNv(); i < n; i++) {
            gradU[0][i] /= angleCount[i];
            gradU[1][i] /= angleCount[i];
        }
    }

    public void computeGradientBySPR(double[] u, double[][] gradU) {
        int len = u.length;
        ParamCheck.checkEqual(len, this.ndof);
        ParamCheck.checkEqual(len, gradU[0].length);
        ParamCheck.checkEqual(len, gradU[1].length);

        double[][] coord = new double[3][2];
        double[][] grad = new double[3][2];
        double[] eleU = new double[3];
        double[][] elementGradient = new double[2][mesh.getNt()];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coord);
            CG.gradient2DT3(coord, grad);
            eleU[0] = u[ele[0]];
            eleU[1] = u[ele[1]];
            eleU[2] = u[ele[2]];
            double[] dx = MVO.col(grad, 0);
            double[] dy = MVO.col(grad, 1);
            elementGradient[0][i] = MVO.dot_product(eleU, dx);
            elementGradient[1][i] = MVO.dot_product(eleU, dy);
        }
        System.out.println("Element Gradient: ");
        System.out.println("dx min: " + MVO.min(elementGradient[0]) +" max: " + MVO.max(elementGradient[0]));
        System.out.println("dy min: " + MVO.min(elementGradient[1]) +" max: " + MVO.max(elementGradient[1]));
        HashSet<Integer>[] patchSet = mesh.SuperconvergencePatches(3);
        int nvIdx = 0;
        for(HashSet<Integer> patch: patchSet) {
            System.out.println(patch);
            int pSize = patch.size();
            double[][] A = new double[pSize][];
            double[] bx = new double[pSize];
            double[] by = new double[pSize];
            int pIdx = 0;
            for(int idx : patch) {
                int[] ele = mesh.getElement(idx);
                mesh.getCoordPerElement(ele, coord);
                double[] xy = mesh.getinteriorPointInElement(coord);
                A[pIdx] = new double[]{1.0, xy[0], xy[1]};
                bx[pIdx] = elementGradient[0][idx];
                by[pIdx] = elementGradient[1][idx];
                pIdx++;
            }
            double[][] AT = MVO.transpose(A);
            double[][] ATA = MVO.matmul(AT, A);
            double[] dx = new double[3];
            double[] dy = new double[3];
            double[] ATBx = MVO.matmul(AT, bx);
            double[] ATBy = MVO.matmul(AT, by);
            DirectDSolver solver = new DirectDSolver(ATA);
            solver.CholeskyFactorize();
            solver.CholeskySolve(ATBx, dx);
            solver.CholeskySolve(ATBy, dy);
            double[] xy = mesh.getNode(nvIdx);
            double[] P = new double[]{1.0, xy[0], xy[1]};
            gradU[0][nvIdx] = MVO.dot_product(P, dx);
            gradU[1][nvIdx] = MVO.dot_product(P, dy);
            nvIdx++;
        }
    }

    public double computeL2Norm(double[] u) {
        double result = 0.0;
        double[] v = new double[4];
        double[] uu = new double[4];
        double[] dx = new double[4];
        double[] dy = new double[4];
        double xi, et, detJ, wxarea, u0;
        for (int k = 0; k < mesh.getNt(); k++) {
            int[] ele = mesh.getElement(k);
            mesh.getCoordPerElement(ele, coordRegion);
            MVO.getSlice(u, ele, uu);
            for (int i = 0; i < 2; i++) {
                xi = GP2[i];
                for (int j = 0; j < 2; j++) {
                    et = GP2[j];
                    detJ = CG.ISOMap2D(CG::Q4ShapeFunction, MVO.col(coordRegion, 0), MVO.col(coordRegion, 1), xi, et, v,
                            dx, dy);
                    u0 = MVO.dot_product(uu, v);
                    wxarea = GW2[i] * GW2[j] * detJ;
                    result += u0 * u0 * wxarea;
                }
            }
        }
        return Math.sqrt(result);
    }

    public double computeH1SemiNorm(double[] u) {
        double result = 0.0;
        double[] v = new double[4];
        double[] uu = new double[4];
        double[] dx = new double[4];
        double[] dy = new double[4];
        double xi, et, detJ, wxarea, dxU, dyU;
        for (int k = 0; k < mesh.getNt(); k++) {
            int[] ele = mesh.getElement(k);
            mesh.getCoordPerElement(ele, coordRegion);
            MVO.getSlice(u, ele, uu);
            for (int i = 0; i < 2; i++) {
                xi = GP2[i];
                for (int j = 0; j < 2; j++) {
                    et = GP2[j];
                    detJ = CG.ISOMap2D(CG::Q4ShapeFunction, MVO.col(coordRegion, 0), MVO.col(coordRegion, 1), xi, et, v,
                            dx, dy);
                    dxU = MVO.dot_product(uu, dx);
                    dyU = MVO.dot_product(uu, dy);
                    wxarea = GW2[i] * GW2[j] * detJ;
                    result += (dxU * dxU + dyU * dyU) * wxarea;
                }
            }
        }
        return Math.sqrt(result);
    }

    public double[] integrate(ScalarFunc[] sfun) {
        int nf = sfun.length;
        double[] result = new double[nf];
        double[][] coef = new double[nf][3];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            int label = mesh.getElementLabel(i);
            getCoefFromFunc(sfun, coordRegion, label, null, coef);

            double area = CG.area2DT3(coordRegion);
            for (int j = 0; j < nf; j++) {
                result[j] += MVO.sum(coef[j]) * area / 3.0;
            }

        }
        return result;
    }

    public double integrate(double[] sol) {
        double result = 0.0;
        double[] coef = new double[3];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            double area = CG.area2DT3(coordRegion);
            coef[0] = sol[ele[0]];
            coef[1] = sol[ele[1]];
            coef[2] = sol[ele[2]];
            result += MVO.sum(coef) * area / 3.0;
        }
        return result;
    }

    public double[] integrate(double[][] sol) {
        int nf = sol.length;
        double[] result = new double[nf];
        double[][] coef = new double[nf][3];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            double area = CG.area2DT3(coordRegion);
            for (int j = 0; j < nf; j++) {
                coef[j][0] = sol[j][ele[0]];
                coef[j][1] = sol[j][ele[1]];
                coef[j][2] = sol[j][ele[2]];
                result[j] += MVO.sum(coef[j]) * area / 3.0;
            }
        }

        return result;
    }

    public void map(double[] source, FES fs, double[] destination) {
        int nv1 = mesh.getNv();
        int nv2 = fs.getMesh().getNv();
        Mesh mesh2 = fs.getMesh();
        ParamCheck.checkEqual(nv1, source.length);
        ParamCheck.checkEqual(nv2, destination.length);
        for (int i = 0; i < nv2; i++) {
            double[] xy = mesh2.getNode(i);
            double interp = getValue(source, xy);
            if (interp != 1.0e30) {
                destination[i] = interp;
            }
        }
    }
}
